//TODO operator=, compilation and testing of all this stuff

#include <inttypes.h>
#include <string>
#include <vector>
#include <cstdlib>
#include <NTL/ZZ.h>
#include "main.h"
#include <sstream>
#include "exceptions.cpp"
#include <sys/time.h>



namespace utl
{
	 
	static const string n_s = "6277101735386680763835789423176059013767194773182842284081";
	static const string p_s = "6277101735386680763835789423207666416083908700390324961279";
	static const string gx_s = "188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012";
	static const string gy_s = "07192b95ffc8da78631011ed6b24cdd573f977a11e794811";			
	static const string zero = "0";

	static int convertNum(char c)
	{
		switch(c)
		{
			case '0':
				return 0;
			case '1':
				return 1;
			case '2':
				return 2;
			case '3':
				return 3;
			case '4':
				return 4;
			case '5':
				return 5;
			case '6':
				return 6;
			case '7':
				return 7;
			case '8':
				return 8;
			case '9':
				return 9;
			case 'A':
			case 'a':
				return 10;
			case 'B':
			case 'b':
				return 11;
			case 'C':
			case 'c':
				return 12;
			case 'D':
			case 'd':
				return 13;
			case 'E':
			case 'e':
				return 14;
			case 'F':
			case 'f':
				return 15;
			default:
				throw charex;
		}
	}
} //end utl namespace


using namespace std;
using namespace NTL;

template <typename Ty>
struct NumDigitsTrait
{
	static const int num_digits;
	static const int lower_bound;
};

template <>
struct NumDigitsTrait<uint8_t>
{
	static const int num_digits = 3;
	static const int lower_bound = 2;
};

template <>
struct NumDigitsTrait<uint16_t>
{
	static const int num_digits = 5;
	static const int lower_bound = 6;
};

template <>
struct NumDigitsTrait<uint32_t>
{
	static const int num_digits = 10;
	static const int lower_bound = 4;
};

template <>
struct NumDigitsTrait<uint64_t>
{
	
	static const int num_digits = 21;
	static const int lower_bound = 1;
};

template <typename Ty>
struct UintTraits
{
	static const uint64_t wordcut = ((uint64_t)1 << (sizeof(Ty) * 8)) - 1;
	
	//General implementation, valid for every type
	//multiplication is done resulting in a uint64_t
	//wich is subsequently cut according to wordsize (defined in main.h)
	static Ty mult(Ty op1, Ty op2, Ty *carry)
	{
		uint64_t uncut_result = ((uint64_t) op1) * ((uint64_t) op2);
		Ty result ;
		if(carry != NULL)	
		{
			*carry = (Ty)(uncut_result >> (sizeof(Ty) * 8));
		}
		uint64_t cut = UintTraits<Ty>::wordcut;
		return result = ((Ty) (uncut_result & cut));
	}
	
	
	static Ty add(Ty op1, Ty op2, Ty *carry)
	{
		Ty result = op1 + op2;
		if(carry != NULL)
			*carry = result < op1;
		return result;
	}
	
	//UNSAFE
	static Ty sub(Ty a, Ty b, Ty *remainder)
	{
		Ty result = a - b;
		(*remainder) -= (result > a);
		return result;
	}

};

// Specialization Trait for uint64_t, 
// whose multiplication exceed most common processors wordsize.
template <>
struct UintTraits<uint64_t>
{
	
	static const uint64_t wordcut = 0xffffffffffffffff;

	// The multiplication needs a substantially different algorithm
	static uint64_t mult(uint64_t op1, uint64_t op2, uint64_t *carry)
	{
		uint64_t result = op1 * op2;
		if(carry != NULL)
		{
			if((op1!= 0) && (result / op1) < op2)
			{
				uint64_t hi_a, lo_a, hi_b, lo_b, cut32, c1;
				cut32 = ((uint64_t)1 << 32) - 1;
				hi_a =  op1 >> 32;
				lo_a =  op1 & cut32;
				hi_b =  op2 >> 32;
				lo_b =  op2 & cut32;
				c1 = (lo_a * lo_b) >> 32;
				c1 += (hi_a * lo_b) & cut32;
				c1 += (lo_a * hi_b) & cut32;
				*carry =  hi_a * hi_b;
				*carry += (hi_a * lo_b) >> 32;
				*carry += (lo_a * hi_b) >> 32;
				*carry += c1 >> 32;
			}
			else
				*carry = 0;
		}
		return result;
	}
	
	//Same as general case, but I have to avoid a compilation error
	static uint64_t add(uint64_t op1, uint64_t op2, uint64_t *carry)
	{
		uint64_t result = op1 + op2;
		if(carry != NULL)
			*carry = result < op1;
		return result;
	}
	
	//UNSAFE
	static uint64_t sub(uint64_t a, uint64_t b, uint64_t *remainder)
	{
		uint64_t result = a - b;
		(*remainder) -= (result > a);
		return result;
	}	

};




template <typename Ty>
class GigInt 
{
	public:
		vector<Ty> chunks;
	
	public:
		// Constructors
		GigInt() 
		{
			chunks.push_back(0);
		}
		
		GigInt(const GigInt<Ty> &that):
			chunks(that.chunks) {}	
	
		GigInt(vector<Ty> &that):
			chunks(that) {}	

		//RISKY. Works only if the two typenames have
		//a dimension multiple of each other
		template <typename Ty2>
		GigInt(vector<Ty2> &that)
		{
			int dim_this = sizeof(Ty);
			int dim_that = sizeof(Ty2);
			int ratio;
			if(dim_this > dim_that)
			{
				ratio = dim_this / dim_that;
				for(int i = 0; i < that.size(); i += ratio)
				{
					Ty chunk = 0;
					for(int j= 0; j < ratio; ++j)
					{
						if((i + j) < that.size())
						{
							chunk = chunk | ((Ty)that[i+j] << (dim_that * 8 * j));
						}
					}
					this->chunks.push_back(chunk);
				}
			}
			else
			{
				ratio = dim_that / dim_this;
				for(int i = 0; i < that.size(); ++i)
				{
					Ty2 data = that[i];
					for(int j = 0; j < ratio; ++j)
					{
						Ty chunk = data & UintTraits<Ty>::wordcut;
						data = data >> dim_this * 8;
						this->chunks.push_back(chunk);
					}
				}
				
			}
		}
		
		GigInt(const char *c_str)
		{
			string str(c_str);
			GigInt<Ty> temp(str);
			this->chunks = temp.chunks;
		}

		GigInt(const ZZ &bignum)
		{
			if(bignum == 0)
			{
				this->chunks.push_back(0);
				return;
			}
			ostringstream o;
		   o << (UintTraits<Ty>::wordcut);
			ZZ zz_cut = to_ZZ(o.str().c_str());
			Ty chunk;
			unsigned long temp2;
			ZZ temp;
			ZZ bignum2 = bignum;
			while(bignum2 > 0)
			{
				temp = bignum2 & zz_cut;
				temp2 = to_ulong(temp);
				chunk = (Ty) temp2;
				this->chunks.push_back(chunk);
				bignum2 = bignum2 >> (sizeof(Ty) * 8);
			}
		}
			
		GigInt<Ty>(string num)
		{
			bool hex = (num.find_first_of("xabcdef") 
				!= string::npos);
			if(hex)
			{
				Ty chunk;
				vector<Ty> results;
				string temp;
				for(int i = num.size() - 1; i >= 0; --i)
					temp.push_back(num[i]);
				int i = 0;
				int finish = num.size();
				if(num[0] == '0' && num[1] == 'x')
					finish = num.size() - 2;
				int step = (sizeof(Ty) * 2);
				for(; i < finish; i += step)
				{
					chunk = 0;
					for(int j = 0; j < step && (i + j) < finish; ++j)
					{
						chunk += (((Ty)utl::convertNum(temp[i+j])) << (4 * j));
					}
					results.push_back(chunk);
				}
				this->chunks = results;
			}
			else
			{
				ZZ temp;
				temp = to_ZZ(num.c_str());
				GigInt<Ty> temp2(temp);
				this->chunks = temp2.chunks;
			}
		}
	

		int size() const
		{
			int valid_chunks = this->chunks.size();
			while(valid_chunks > 0 && (this->chunks[valid_chunks-1] == 0))
				--valid_chunks;
			return valid_chunks;
		}

		int bitSize()const
		{
			return this->size() * sizeof(Ty) * 8;
		}


		Ty intPow (Ty base, int exp, int *err)
		{
			if(exp == 0) return base;
			Ty result, old_result = 1;
			for(int i=0;i<exp;++i)
			{
				result = old_result * 10;
				if (result%10 !=0)
				//if (result/old_result != 10)
				//if(old_result > result)
				{
					*err = 1;	
					return 0;
				}
				old_result = result;
			}
			Ty res, carry;
			res = UintTraits<Ty>::mult(result, base, &carry);
			if(carry)
			{
				*err = 1;
				return 0;
			}
			return res;
		}
		
		void dumpDebug() const
		{
			
			for(int i=chunks.size()-1;i>=0;--i)
			{
				if(chunks[i])
				{
					cerr << (uint64_t)chunks[i];
					cerr << " ";
				}
				else 
				{
					int num_digits;
					if((i < chunks.size()-1) && 
							(chunks[i+1] & 1 > NumDigitsTrait<Ty>::lower_bound))
						num_digits = NumDigitsTrait<Ty>::num_digits - 1;
					else
						num_digits = NumDigitsTrait<Ty>::num_digits;
					for(int i = 0; i < num_digits; ++i)
						cerr << "0";
					cerr << " ";
				}
			}	
			cerr << endl;
		}

		string to_str() const
		{
			ZZ temp;
			temp = GigInt::GtoZZ(*this);
			stringstream result;
			result << temp;
			return result.str();
		}

		void dump() const 
		{
			cerr << (*this).to_str()  << endl;
		}
		
		
		GigInt operator=(const GigInt &that)
		{
			this->chunks = that.chunks;
			return *this;
		}		

		//Extension
		static GigInt<Ty> extendTo(const GigInt<Ty> &a, int new_size)
		{
			vector<Ty> results(a.chunks);
			int old_size = results.size();
			for(int i = 0; i < (new_size - old_size); ++i)
				results.push_back((Ty)0);
			GigInt<Ty> res(results);
			return res;
		}	

		//Comparation
		bool operator==(const GigInt &that) const
		{
			GigInt<Ty> op1, op2;
			GigInt<Ty>::makeEqual(*this, that, &op1, &op2);
			for(int i = 0; i < op1.chunks.size(); ++i)
			{
				if(op1.chunks[i] != op2.chunks[i])
					return false;
			}
			return true;
		}
		bool operator!=(const GigInt &that) const
		{ return !(*this == that); }


		bool operator<(const GigInt &that) const
		{
			ZZ op1, op2;
			op1 = GigInt<Ty>::GtoZZ(*this);
			op2 = GigInt<Ty>::GtoZZ(that);
			return (op1 < op2);
		}

		bool operator>(const GigInt &that) const
		{
			ZZ op1, op2;
			op1 = GigInt<Ty>::GtoZZ(*this);
			op2 = GigInt<Ty>::GtoZZ(that);
			return (op1 > op2);
		}

		bool operator>=(const GigInt &that) const
		{
			return !(*this < that);
		}

		bool operator<=(const GigInt & that) const
		{
			return !(*this > that);
		}
		

		// Addition
		GigInt operator+(const GigInt<Ty> &that) const
		{
			vector<Ty> results;
			vector<Ty> op1(this->chunks);
			vector<Ty> op2(that.chunks);
			Ty res_chunk;
			Ty overflow=0;
			int n, n_this, n_that;
			n_this = op1.size();
			n_that = op2.size();
			if(n_this == n_that) //most likely
				n = n_this;
			else
			{
				GigInt<Ty> temp;
				if(n_this > n_that)
				{
					n = n_this;
				   temp = GigInt<Ty>::extendTo(that, n_this);
					op2 = temp.chunks;
				}
				else
				{
					n = n_that;
					temp = GigInt<Ty>::extendTo(*this, n_that);
					op1 = temp.chunks;
				}
			}
			Ty k_in, k_out = 0, overflow2;		
			for(int i = 0; i < n; ++i)
			{
				k_in = k_out;
				res_chunk = UintTraits<Ty>::add(op1[i], op2[i], &overflow);
				res_chunk = UintTraits<Ty>::add(res_chunk, k_in, &overflow2);
				k_out = overflow + overflow2;	
				results.push_back(res_chunk);
			}
			
			if(k_out)
			{
				results.push_back(k_out);
			}
			
			GigInt<Ty> res(results);
			return res;
		}
		
		GigInt operator++ ()
		{
			string s_one = "1";
			GigInt<Ty> one(s_one);
			*this = *this + one;
			return *this;
		}

		
		inline GigInt operator*(const GigInt<Ty> &that) const
		{
			return GigInt<Ty>::karatsubaMul(*this,that);
		}


		

		//Subtraction
		GigInt operator-(const GigInt<Ty> &that) const
		{
			
			if((*this) < that)
				throw subex;
			GigInt<Ty> op1, op2;
			GigInt<Ty>::makeEqual(*this, that, &op1, &op2);
			vector<Ty> results(op1.size(), 0);
			Ty remainder = 0;
			int size = op1.size();
			results[size - 1] = UintTraits<Ty>::sub(op1.chunks[size -1], 
					op2.chunks[size -1], &remainder);
			for(int i = op1.size() - 2; i >= 0; --i)
			{
						
				Ty oldRemainder = results[i + 1];
				results[i] = UintTraits<Ty>::sub(op1.chunks[i], op2.chunks[i], 
						&results[i + 1]);
				if(results[i + 1] == (Ty)-1 && results[i + 1] != oldRemainder)
				{
					bool exit = false;
					for(int j = i + 2; !exit && (j < size); ++j)
					{
						if(results[j] != 0)
							exit = true;
						--results[j];
					}
				}
			}
			GigInt<Ty> result(results);
			return result;
		}

		
		// Multiplications
		// WARNING
		// We are assuming that the two operands have same length
		static GigInt schoolbookMul(const GigInt<Ty> &a, const GigInt<Ty> &b)
		{
			vector<Ty> results(a.chunks.size()*2, 0);
			int n = a.chunks.size();
			Ty k, u, v;
			for (int j = 0;  j < n; ++j)
			{
				k = 0;
				for (int i = 0; i < n; ++i)
				{
					v = UintTraits<Ty>::mult(a.chunks[i], b.chunks[j], &u);
					v += k;
					u += (v < k);
					v += results[i+j];
					u += (v < results[i+j]);
					results[i+j] = v;
					k = u;
				}
				results[j+n] = k;
			}
			return GigInt(results);
		}

		//TODO Check if operands have different sizes
		static GigInt schoolbookMulFaulted(const GigInt<Ty> &a, 
				const GigInt<Ty> &b, int ind)
		{
			int i_fault = ind, j_fault = ind;
			GigInt<Ty> temp_a = a;
			Ty random_val = temp_a.chunks[ind] + (Ty)rand();

			vector<Ty> results(a.chunks.size()*2,0);
			int n = a.chunks.size();
			Ty k, u, v;
			for (int j = 0;  j < n; ++j)
			{
				k = 0;
				for (int i = 0; i < n; ++i)
				{	
					if(i == i_fault && j == j_fault)
						v = UintTraits<Ty>::mult(random_val, b.chunks[j], &u);
					else
						v = UintTraits<Ty>::mult(a.chunks[i], b.chunks[j], &u);
					v += k;
					u += (v < k);
					v += results[i+j];
					u += (v < results[i+j]);
					results[i+j] = v;
					k = u;
				}
				results[j+n] = k;
			}
			return GigInt(results);
		}

		//Here we make the operands have the same length
		static GigInt karatsubaMul (const GigInt<Ty> &a, const GigInt<Ty> &b)
		{
			GigInt<Ty> op1, op2;
			GigInt<Ty>::makeEqual(a, b, &op1, &op2);
			int MIN_LEN = 320;	
			int len = op1.bitSize(); //min(a.bitSize(),b.bitSize());
			//int chunks = max(a.chunks.size(),b.chunks.size());
			if (len <= MIN_LEN)	
				return schoolbookMul(op1, op2);

			GigInt<Ty> temp_a = op1;
			GigInt<Ty> temp_b = op2;

			/** Not necessary if both num have the same length
			if(a.chunks.size() < chunks)
				temp_a = GigInt<Ty>::extendTo(temp_a,temp_b.chunks.size());
			else if (b.chunks.size() < chunks)
				temp_b = GigInt<Ty>::extendTo(temp_b,temp_a.chunks.size());
			**/

			if(temp_a.chunks.size()%2!=0)
			{
				temp_a = GigInt<Ty>::extendTo(temp_a,temp_a.chunks.size()+1);
				temp_b = GigInt<Ty>::extendTo(temp_b,temp_b.chunks.size()+1);
			}
			
			vector<Ty> lo_a(temp_a.chunks);
			lo_a.resize(lo_a.size()/2);
			vector<Ty> lo_b(temp_b.chunks);
			lo_b.resize(lo_b.size()/2);
	
			vector<Ty> hi_a(temp_a.chunks);
			reverse(hi_a.begin(), hi_a.end());
			hi_a.resize(hi_a.size()/2);
			reverse(hi_a.begin(), hi_a.end());
			vector<Ty> hi_b(temp_b.chunks);
			reverse(hi_b.begin(), hi_b.end());
			hi_b.resize(hi_b.size()/2);
			reverse(hi_b.begin(),hi_b.end());

				
			GigInt<Ty> g_lo_a = GigInt(lo_a);
			GigInt<Ty> g_hi_a = GigInt(hi_a);
			GigInt<Ty> g_lo_b = GigInt(lo_b);
			GigInt<Ty> g_hi_b = GigInt(hi_b);
			
			GigInt<Ty> g_sum_a = g_lo_a + g_hi_a;	
			GigInt<Ty> g_sum_b = g_lo_b + g_hi_b;	

			GigInt<Ty> g_hi_hi = karatsubaMul(g_hi_a,g_hi_b);
			GigInt<Ty> g_lo_lo = karatsubaMul(g_lo_a,g_lo_b);
			GigInt<Ty> g_sum_a_b = karatsubaMul(g_sum_a,g_sum_b);

			//cout << "g_hi_hi "; g_hi_hi.dump(); 
			//cout << "g_lo_lo "; g_lo_lo.dump(); 
			//cout << "g_sum_a_b "; g_sum_a_b.dump();
			GigInt<Ty> g_z1 = g_sum_a_b - g_hi_hi - g_lo_lo;	
			//cout << "g_z1 "; g_z1.dump();

			return (shiftLeft(g_hi_hi,(temp_a.chunks.size())) + 
					shiftLeft(g_z1,(temp_a.chunks.size()/2)) + g_lo_lo);
		}

		
		GigInt operator* (const Ty &that)
		{
			vector<Ty> results;
			vector<Ty> res_round;
			Ty carry, temp;
			for(int i = 0; i < this->chunks.size(); ++i)
			{
				temp = UintTraits<Ty>::mult(this->chunks[i], that, &carry);
				res_round.push_back(temp);
				if(carry)
					res_round.push_back(carry);
				for(int k = 0; k < res_round.size(); ++k)
				{
					if(results.size() > k+i)
						results[k+i] += res_round[k];
					else
						results.push_back(res_round[k]);
				}
			}
			return GigInt<Ty>(results);
		}

		//Division
		//Yet to completely implement
		//It's only a wrapper to a ZZ division
		GigInt operator/(const GigInt<Ty> &b) const
		{
			ZZ op1, op2;
			op1 = GigInt::GtoZZ(*this);
			op2 = GigInt::GtoZZ(b);
			GigInt<Ty> result(op1 / op2);
			return result;
		}
		//Same as operator/
		GigInt operator%(const GigInt<Ty> &that) const
		{
			ZZ op1, op2;
			op1 = GigInt::GtoZZ(*this);
			op2 = GigInt::GtoZZ(that);
			GigInt<Ty> result(op1 % op2);
			return result;
		}
		

		static ZZ GtoZZ(const GigInt &gig)
		{
			ZZ result = to_ZZ(0);
			ZZ temp;
			unsigned long exponent = sizeof(Ty) * 8;
			for(int i = 0; i < gig.chunks.size(); ++i)
			{
				temp = to_ZZ((unsigned long) gig.chunks[i]);
				temp = temp << (i * exponent);
				result += temp;
			}
			return result;
		}


		//This functions has the effect to multiply a GigInt 
		//by a 2^wordsize*words factor
		static GigInt shiftLeft(const GigInt<Ty> &a, int words)
		{
			vector<Ty> results;
			for(int i = 0; i < words; ++i)
				results.push_back((Ty)0);
			for(int j = 0; j < a.chunks.size(); ++j)
				results.push_back(a.chunks[j]);
			return GigInt<Ty>(results);
		}

		static void makeEqual(const GigInt<Ty> &a, const GigInt<Ty> &b,
			  	GigInt<Ty> *op1, GigInt<Ty> *op2)
		{
			if(op1 == NULL || op2 == NULL)
				return;
			if(a.size() > b.size())
			{
				*op1 = a;
				*op2 = GigInt<Ty>::extendTo(b, a.size());
			}
			else
			{
				if(a.size() < b.size())
				{
					*op1 = GigInt<Ty>::extendTo(a, b.size());
					*op2 = b;
				}
				else
				{
					*op1 = a;
					*op2 = b;
				}
			}
		}



void randGigInt(string min, string max)
{
	GigInt<Ty> maxVal(max);
	GigInt<Ty> minVal(min);
	
	ZZ upperBound = to_ZZ(max.c_str());
	ZZ result = to_ZZ(0);
	struct timeval t;
	long usec;
	gettimeofday(&t, NULL);
	usec = t.tv_usec;
	SetSeed(to_ZZ(usec));	
	while(result==0)
	{
		result = RandomBnd (upperBound);
	}

	GigInt<Ty> res(result); 
	this->chunks = res.chunks; 
}


GigInt mod(string m)
{
	GigInt<Ty> mod(m);
	
	while(*this > mod)
		*this = (*this) - mod;
	return  *this;
}




void  mod192(){
	GigInt<uint64_t> A(this->chunks);
	if(A.size()<6)
		A = GigInt<uint64_t>::extendTo(A,6);
	GigInt<uint64_t> T;
	T.chunks.pop_back();
	T.chunks.push_back(A.chunks[0]);
	T.chunks.push_back(A.chunks[1]);
	T.chunks.push_back(A.chunks[2]);

	GigInt<uint64_t> S1;
	S1.chunks.pop_back();
	S1.chunks.push_back(A.chunks[3]);
	S1.chunks.push_back(A.chunks[3]);

	GigInt<uint64_t> S2;
	S2.chunks.push_back(A.chunks[4]);
	S2.chunks.push_back(A.chunks[4]);

	GigInt<uint64_t> S3;
	S3.chunks.pop_back();
	S3.chunks.push_back(A.chunks[5]);
	S3.chunks.push_back(A.chunks[5]);
	S3.chunks.push_back(A.chunks[5]);

	GigInt<uint64_t> B = T + S1 + S2 + S3;
	GigInt<Ty> result(B.chunks);
	GigInt<Ty> p(utl::p_s);
	if(result >= p)
	{
		ZZ num_z = GigInt<Ty>::GtoZZ(result);
		ZZ p_z = GigInt<Ty>::GtoZZ(p);
		ZZ_p::init(p_z);
		ZZ_p num_p = to_ZZ_p(num_z);
		GigInt<Ty> temp(NTL::rep(num_p));
		this->chunks = temp.chunks;
		return;
	}
	this->chunks = result.chunks;
}

};

